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Figure  Captions 

Figure  1.  Flow  chart  showing  the  major  elements  of  the  ionospheric  chemistry 
model. 

Figure  2.  Electron  density  profiles  from  0-60  minutes  using  precipitating 
electron  flux  from  the  PIIE  rocket  experiment;  local  chemistry  only  (no  0+ 
diffusion). 

Figure  3.  Same  as  Figure  2  but  including  0+  diffusion. 

Figure  4a-d.  Same  as  Figure  2  (0-10  minutes)  but  with  high  altitude  boundary 
0+  flux  set  at  1  x  lO1^  m'2  s"1  (down),  1  x  10^  (down),  1  x  101®  (up),  and  5  x 
10*2  (Up),  respectively. 

Figure  5.  Schematic  showing  the  contribution  to  production  by  flux=<j>i  at 
energy=E5. 

Figure  6a-f.  Altitude  profiles  of  production  rates  for  02+,  0+,  and  N2+ 
contributed  by  <J>1=1  x  10^  (cm‘2  sr‘l  eV'l)  at  E=55,  180,  320,  550,  750,  and 
1300  eV. 

Figure  7.  Differential  electron  flux  measured  by  the  AFGL  PIIE  rocket 
experiment  within  an  auroral  arc  on  March  15,  1985. 

Figure  8a.  Maxwellian  flux  with  200  eV  characteristic  energy. 

Figure  8b.  Production  rate  profiles  calculated  from  the  tables  for  the  flux  in 
Figure  8a. 

Figure  9.  PIIE  observations  of  electron  density  fluctuations,  total  electron 

density,  electron  temperature,  and  energy  flux  of  precipitating  electrons  vs. 
time  (seconds  after  launch),  altitude,  and  UT.  Vertical  line  indicates  time  of 
the  model  calculation  inside  the  arc. 

Figure  10.  Model  production  rates  of  0+(4S),  N+  and  N2+.  Also  shown  is  the 
plasma  heating  rate  normalized  to  the  ambient  electron  density.  The  curves 
labelled  (a)  and  (b)  refer  to  the  electron  density  profiles,  shown  in  Figure  12, 

that  are  used  in  the  flux  transport  code.  Production  rates  corresponding  to 
each  curve  are  nearly  identical. 

Figure  11.  Ten  point  averages  of  PIIE  electron  differential  flux  measurements 

(downgoing  electrons).  Sample  from  282.3  to  288.1  s  is  used  to  model  inside  the 
arc.  Other  sample  is  used  to  establish  an  initial  electron  density  profile  for  the 
ionospheric  chemistry  model. 

Figure  12.  Electron  density  and  Te  results  from  the  chemistry  model  using  the 

rates  shown  in  Figure  10.  Densities  shown  correspond  to  rates  calculated  for 
curve  (a). 

Figure  13.  Hall  and  Pedersen  electric  conductivities  for  the  5  minute  profile  in 
Figure  12. 
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Figure  14.  Calculated  PIIE  electron  density  profiles  corresponding  to  flux 
durations  of  1-30  minutes.  Initial  0+  density  is  104  cm '3  at  all  altitudes. 

Figure  15.  Same  as  Figure  14  except  starting  from  a  lower  initial  electron 
density  (10^  cm‘3). 

Figure  16.  Decay  of  electron  density  profile  after  flux  is  turned  off. 

Figure  17.  Development  of  electron  density  profile  using  electron  flux  outside 
arc . 

Figure  18.  Same  as  Figure  17  except  starting  from  a  lower  initial  electron 
density  (10^  cm'3). 

Figure  19.  Calculated  electron  density  profiles  using  flux  inside  arc  and  a  cold 
neutral  atmosphere  (Tn=750  K). 

Figure  20.  Same  as  Figure  19  except  with  a  hotter  neutral  atmosphere 
(Tn=1200K). 
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1.  Introduction 


A  realistic  model  of  the  ionosphere  is  an  extremely  complicated  physical 
and  mathematical  problem.  Excitation  and  ionization  by  precipitating 
particles  or  solar  radiation  influence  chemical  reactions  at  every  altitude. 
Diffusion  along  field  lines  and  hr  vontai  transport  determine  the  final 
electron  and  ion  densities  at  any  point.  These  local,  unsteady  processes  are 

superimposed  on  global  and  secular  variations  having  diurnal,  seasonal,  and 
solar-cycle  dependences. 

In  the  auroral  zone  and  polar  cap  it  is  clear  that  precipitating  electrons  and 
protons  have  a  profound  effect  on  the  ionospheric  plasma  density.  However, 
electric  fields  may  convect  F-region  plasma  hundreds  or  thousands  of 
kilometers  from  its  source;  hence,  optical  emissions  and  density  fluctuations 
associated  with  large  scale  plasma  structures  can  also  occur  some  distance 
away  from  direct  sources  of  ionization.  Although  satellites,  rockets,  radars, 

and  photometers  have  probed  ionospheric  structures  in  detail,  the  origin  and 
fate  of  plasma  in  the  polar  cap  remain  difficult  problems.  A  reliable  model  of 
ionospheric  response  to  ionizing  sources  would  help  to  unravel  some 
difficulties  by  defining  the  portion  of  ionization  produced  locally. 

Many  researchers  have  modeled  the  ionosphere  with  varying  degrees  of 
complexity  depending  primarily  on  the  application  (Jones  and  Rees.  1973; 

Roble  and  Rees.  1977;  Schunk  et  al..  1986).  Global  descriptions  under  different 
geophysical  conditions  rely  heavily  on  large  scale  empirical  studies;  physical 
models  based  on  the  governing  equations  have  been  limited  by  the  availability 
of  key  parameters  and  computational  difficulties.  The  ability  to  test  and 
refine  a  model  has  grown  with  the  aid  of  coordinated  experiments, 
comprehensive  data  sets,  and  large,  high-speed  computers. 
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The  model  developed  by  D.  Strickland  and  colleagues  has  two  distinct  parts, 
flux  transport  and  ionospheric  chemistry.  Flux  transport  begins  with  a 
specified  differential  electron  flux  at  some  upper  altitude  boundary  and 
calculates  electron  fluxes  throughout  the  ionospheric  depth.  This  is 
accomplished  by  a  fast  and  accurate  solution  of  the  flux  transport  equation 
that  treats  both  elastic  and  inelastic  scattering  by  the  principal  neutral 
species,  N2,  O2,  and  O.  A  table  of  cross  sections  for  inelastic  scattering  and 

neutral  densities  describe  the  details  of  each  interaction.  An  eigenvalue 
representation  solves  for  differential  fluxes  over  a  range  of  altitudes  and 
energies  at  20  pitch  angles  from  0  to  180  degrees.  Production  rates  for  selected 
ionic  and  neutral  species  are  found  from  the  calculated  fluxes  and  tabulated 
cross  sections  for  ionization  and  various  excitations.  The  original  code  with 
electron  fluxes  has  been  extended  by  Strickland  to  include  both  precipitating 
protons  and  solar  EUV  input;  results  in  this  report  come  from  earlier  versions 
with  electron  fluxes  only. 

Production  rates  from  the  transport  part,  whatever  the  source  of  ionization, 
become  inputs  to  the  ionospheric  chemistry  model.  Our  attention  has  centered 
on  modifying  the  chemistry  programs  to  handle  production  generated  by  soft 
(1  keV)  fluxes  of  electrons  typically  observed  over  polar  cap  auroral  arcs.  This 
interest  required  some  changes  inherent  in  F-region  physics,  chiefly  by  the 
integration  of  an  existing  0+  diffusion  code  (written  by  Rob  Danniel)  with 
local  chemistry  and  the  addition  of  new  electron  and  ion  temperature 
calculations.  The  following  sections  describe  these  codes  with  some 
applications  and  discuss  remaining  problems  and  goals. 

The  major  scientific  application  of  this  study  can  be  found  in  Weber  et  al. 
(1988),  which  analyzes  results  from  the  PIIE  rocket  experiment  conducted  in 
Greenland  on  March  15,  1985.  To  complement  dc  electric  field  measurements 
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taken  during  the  flight.  Hall  and  Pedersen  electric  conductivity  calculations 
were  added  to  the  model  based  on  standard  expressions  involving  calculated 
electron  densities  and  temperatures.  Combining  this  information  with  the 
measured  electric  fields,  a  separate  program  determined  field-aligned  currents 
from  the  continuity  equation  and  tensor  conductivity. 

In  a  separate  effort  in  support  of  the  AFGL  airborne  ionospheric 
observatory,  J.  Fleischman  participated  in  several  missions  and  became 
involved  in  improving  the  collection  system  of  navigational  data  needed  for 
mission  analysis.  A  Zenith  Z-248  laptop  conputer  was  programmed  for  storage 
of  INS  data  to  back  up  the  existing  tape  system  and  for  a  running  presentation 
of  current  navigational  data  and  related  navigational  problems.  A  manual 
describing  the  system  and  its  operation  was  written  by  J.  Fleischman  as  a 
technical  report. 

The  bulk  of  this  report  centers  on  mathematical  details  and  physical 
assumptions  in  the  model  code  because  the  model  is  essentially  a  tool  which 
should  be  applied  thoughtfully  in  light  of  its  limitations.  In  any  case,  a 
knowledgable  user  will  eventually  become  familiar  with  these  details  and  will 
adjust  the  internal  workings  of  the  model  and  the  peripheral  procedures  to 
suit  his  own  purposes  and  tastes. 
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2.  Details  of  the  Model 


A.  Overview 

Designed  for  E-region  applications,  the  original  code  solved  local  chemical  rate 
equations  for  the  species  listed  in  Table  1.  Later,  0+  diffusion  was  added  to  extend  the 
altitude  range  into  the  F-region.  We  made  structural  changes  in  the  programs  to 

allow  the  chemical  rate  and  diffusion  coefficients  to  reflect  density  information  from 
each  part.  In  addition,  instead  of  assuming  a  Fixed  electron  temperature  (Te)  profile, 

we  developed  a  code  to  solve  the  electron  heat  conduction  equation  based  on  the 
computed  electron  density.  The  eiectron  plasma  heating  rate  is  supplied  by  the  same 
flux  transport  code  that  calculates  chemical  production  rates. 

The  flow  chart  in  Figure  1  illustrates  the  procedure  followed  by  the  model  as  it 
now  stands.  The  lower  half  pertains  to  ionospheric  chemistry  and  thermal  energy 

balance.  Inputs  are  production  rate  (cm"3  s"^).  Starting  with  initial  density  profiles 

for  all  the  species  and  an  initial  Te  profile,  the  code  integrates  the  quantities 
(densities  and  Te)  until  some  desired  time  interval  has  elapsed  or  equilibrium  is 
reached.  Since  diffusion  and  heat  conduction  are  represented  by  partial  differential 

equations  in  time  and  altitude,  appropriate  boundary  conditions  are  also  required. 

Because  Te  responds  very  quickly  to  changes  in  the  heating  rate,  time  steps  must 

be  very  short  at  First  to  handle  potentially  large  temperature  variations.  During  the 
First  second  after  t=0,  a  cycle  through  the  basic  program  unit  (THERM,  CHEM,  and  0+ 
DIFFUS)  is  performed  every  .1  second.  From  t=I  to  120  s  the  update,  dt,  is  lengthened 
to  1  s;  after  120  s,  dt=I0  s.  These  times  are  appropriate  for  the  types  of  conditions  we 
have  encountered,  but  since  the  internal  logic  of  each  unit  is  the  same,  regardless 
the  time  increment  dt,  experimentation  is  not  difficult.  In  fact,  each  program  unit 
has  its  own  internal  time  step  (calculated  in  CHEM  and  0+  DIFFUS,  set  in  THERM),  so 
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that  the  overall  update  time  dt  only  affects  the  overhead  incurred  in  recomputing 
various  coefficients. 

At  the  start  of  an  arbitrary  time  interval,  tn,  THERM  calci'ates  the  thermal 
conduction  coefficient  for  electrons  from  the  current  values  of  electron  density,  nc, 
and  Te.  It  then  solves  for  Te  at  a  time  ln+  1  =tn  +  8t  using  coefficients  which  are 
assumed  constant  during  St.  The  new  Te  resets  temperature  dependent  chemical  rate 

coefficients  for  density  calculations  in  CHEM.  When  all  densities  except  0+  are  found 
at  time  tn+i,  the  new  chemical  rate  (from  CHEM)  and  Te  dependent  diffusion 
coefficient  for  0+  go  to  0+  D1FFUS,  which  finds  the  0+  density  at  tn+j.  The  resulting 
total  electron  density  (nc=nj)  and  Te  at  t-tn+ 1  now  become  inputs  to  THERM  for  the 
next  cycle  beginning  at  tn+j.  By  requiring  that  the  overall  update  time,  dt,  reflects 
the  fastest  changing  quantity,  in  this  case  Te,  we  can  minimize  errors  introduced  by 
solutions  with  constant  coefficients  over  discrete  time  intervals.  A  more  elaborate 
treatment  would  entail  representative,  average  values  of  the  coefficients  and  an 
update  time  calculated  to  ensure  a  specified  accuracy. 

At  present  both  the  electron  flux  transport  and  chemistry  codes  are  one 
dimensional  with  dependences  only  along  a  magnetic  field  line.  Given  the 
complexity  of  the  existing  transport  code,  this  most  likely  will  remain  a  practical 

limitation;  however,  if  the  electron  gyroradius  is  small  compared  with  the  horizontal 
extent  of  the  incident  flux,  true  for  most  auroral  arcs,  the  approximation  is 
acceptable.  Chemical  lifetimes,  on  the  other  hand,  are  comparable  to  transport  times 

of  F-region  density  structures  by  convection  electric  fields.  When  this  occurs 
integration  times  must  be  held  small  or  convection  velocities  must  be  assumed  to  lie 
along  a  constant  density  structure  like  an  arc. 

Temporal  variations  pose  a  somewhat  different  problem.  Although  flux  transport 
is  a  steady  state  formulation,  the  time  required  to  establish  production  rates  arc  short 

(Is)  and  effectively  vary  with  the  input  flux.  A  separate  flux  transport  calculation 
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would  be  required  for  each  new  incident  flux,  not  a  practical  undertaking  in  most 
cases.  However,  if  spectral  shapes  are  similar,  scaling  the  production  rates  relative 
to  an  absolute  flux  could  mimic  temporal  flux  variations.  Production  rates  passed  to 

the  chemistry  code  arc  held  constant,  or,  after  operating  for  a  designated  time,  can 
be  switched  off  to  simulate  decay  of  ion  density.  There  is  no  reason  why  direct 
production  by  an  incident  flux  must  be  fixed  in  CHEM;  the  only  requirement  is  that 
variations  be  slow  compared  with  the  time  steps  in  the  integrating  routines. 

Each  basic  code  unit,  CHEM,  0+  DIFFUS,  and  THERM  involves  the  continuity 
equation  for  scalar  quantity,  ion  density  (internal  energy),  of  the  form 

H=V.a  +  q-l 

(2.1) 

where  f  is  ion  (energy)  density  in  ion  (eV)  cm"3,  S  is  a  flux  in  ions  (eV)  in  cm'^  s'* 
and  q  and  1  are  source  and  loss  terms,  respectively,  in  ions  (eV)  cm‘3  s'*.  In  one 
dimension  the  only  spatial  variation  is  altitude,  or  more  accurately  along  a  magnetic 
Held  line  where  dh=ds  sinl,  h=altitude,  s=distance  along  B  (magnetic  field),  and 
I=inclination  of  B.  Near  the  poles  B  is  nearly  vertical  and  dh>ds  which  we  assume 
here. 

B.  Chemistry 

Since  diffusion  is  relatively  unimportant  in  the  lower  ionosphere,  densities  of 
various  species  are  governed  by  a  system  of  local  chemical  rate  equations, 

3N' 

yf  =  Pi  +  rr^n  kmn  Nm  Nn  -  Nj  ^  lm  Nm  -  ai  Nj  ^  Nm 

(2.2) 
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Pi  represents  direct  production  by  any  process,  e.g.,  electron  impact;  kmn  arc  rate 
coefficients  for  reactions  involving  species  m  and  n  that  produce  species  i. 

Similarly,  lm  and  aj  describe  loss  reactions  with  species  m  and  recombination  with 
electrons  (ne=  L  nj)  ,  respectively.  The  rate  coefficients  k,  1,  and  a  are  read  from  an 

input  data  file  when  the  program  starts.  We  have  modified  some  rates  and  added 
reactions  to  the  set  supplied  with  the  original  Strickland  code  based  on  a  table  in 
VallaneeJones  (19741.  The  rates  can  be  changed  directly  in  the  input  file  or  through 
the  Cyber  UPDATE  utility.  Rates  depending  on  the  electron  temperature  arc 
interpolated  from  a  table  computed  with  the  current  Te  profile. 

Subroutine  CHEMEQ  performs  the  integration  of  (2.2)  for  each  species  over  a  time 
step  beginning  at  T  and  ending  at  T+DTUP.  Starting  with  densities  calculated  from 
the  previous  step,  a  preparatory  routine,  DERIV,  computes  the  necessary  coefficients 
and  calculates  densities  for  those  species  designated  as  "fast",  which  adjust  quickly  to 
changes  in  production  and  loss  and  hence  are  close  to  equilibrium  at  all  times. 
Assuming  0Ni/9ti  =  0  for  the  fast  species,  Nj  =  pj/qj.  The  species  list  (Table  1)  is 

arranged  so  that  all  the  "slow"  species  to  be  integrated  appear  first  starting  with 
index  5,  after  which  the  "fast"  species  are  listed.  Unless  otherwise  noted,  species  5 
through  10  (NO+,  O2+,  0+,  N(4S),  and  N(2D))  are  numerically  integrated  in  the 

examples  shown  in  this  report.  Species  1 1  through  20  are  assumed  to  be  in 
equilibrium  in  the  sense  described  above. 

The  numerical  integration  routine  CHEMEQ  employs  a  second  order  predictor- 
corrector  scheme.  Letting  Y  represent  the  density,  N  ,  we  can  express  (2.2)  as 


dY 

dt 


=  C  -  D 
(2.3) 
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where  C  is  production  and  D  =  Y/x  stands  for  the  loss  term  which  is  proportional  the 
independent  variable,  Y.  Using  an  internal  time  step  t,  chosen  to  reflect  the  fastest 
changing  species,  a  finite  difference  representation  of  (2.3)  involves  quantities  at 
the  initial  time,  to,  and  final  time,  tQ+t  . 


Y-Y0 


=  C0  - 


t 


o 


(2.4) 


The  denominator  in  the  LHS  of  (2.4)  emphasizes  that  the  finite  difference  is  a  second 
order  representation  of  the  time  derivative  at  the  midpoint;  accordingly,  the  RHS  is 
an  average  (C  =  CG,  and  x  =  t0). 

From  (2.4)  a  predicted  value  of  Y  after  a  time  interval  t  can  be  found, 


Y o  (2  T0  -  t)  +  2x0tc0 

(2  x0  +  t) 


(2.5) 


With  the  predicted  value  of  Yp  new  average  values  of  x  and  C  can  be  calculated  to 
yield  a  corrected  value  Y  after  another  application  of  (2.5). 


x  = 


(Yo  Ynl  1  , 

[d^+  D^J*2<x0+‘) 


C  =  g  (Co  ■  Cp) 
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[Y0(2x -t)  +  2X  +  C] 

Y  c  = 

(2*  +  t 

(2.6) 

Predicted  values  for  C  and  t  at  t  =  to+t  are  evaluated  by  a  call  to  subroutine  DERIV.  If 
Y  c  fails  a  convergence  criterion  the  time  step  t  is  reduced  and  the  procedure 
repeated. 

Up  to  order  t^  ,  the  integration  method  in  CHEMEQ  is  identical  to  results  obtained 
from  a  standard  modified  Euler  numerical  integration.  Using  the  same  terminology 
as  above,  the  predicted  Euler  value  is  simply  a  linear  extrapolation  from  the  initial 
value. 


Yp  =  Y0  +  Y0' t  *  Y0  + 


t 


(2.7) 


A  corrected  slope  is  obtained  from  an  average  with  the  predicted  value. 


Y'ave-J(Y0'  +  Yp-) 


"" 

L  yo) 

f  Y0)  Y°  + 

co  -  + 

co  -  7" 

l  T°j 

t 

(2.8) 


The  corrected  value  then  becomes, 
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YcEuler  =  y0  + 

V'ave  t 


V0 


2  x, 


2t  + 


t2  ^ 


2x, 


+  Cq  t 


(2x0  -  t) 


2x, 


(2.9) 


Expanding  the  denominator  in  (2.5)  by  a  Taylor  series  we  get. 


2t0  +  t 


(A  J_Y  i 
2xrt 


2t, 


_L  t2 

2t0  +  (2x0)2  + 


so  that. 


(2.10) 


Yp  =  Y0  (2x0  - 1) 


1  -L-  *2  N 

2*0  *  (2x0)2 

2x„ 


+  Co2xnt 


f  j_  _ti_( 

‘  2t0  +  (2t0)2 


2x, 


(2.11) 


or. 


Vo 


2x, 


t2^ 
2t  +  — 


2t, 


Y°  4x03  C0t  (2xp  -  t)  C°2x0 

2xq  2to  2xq 
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=  ycEuler  +  0(t3) 

(2.12) 


C.  Diffusion 

General  discussions  about  vertical  transport  in  the  ionosphere  can  be  found  in 
Rishbeth  and  Garriot  (1969)  and  Banks  and  Kockarts  (1973).  The  starting  point  is  an 
expression  for  the  diffusion  velocity  Wj  of  an  ion  along  a  magnetic  field  line, 


Wd 


=  -Da  sin2! 


1_d_N 
N  dz  + 


HpJ 


(2.13) 


Da  is  the  plasma  diffusion  coefficient,  N  the  ion  density,  z  the  distance  along  a  field 
line,  Tp=(Tj+Te)/2  the  plasma  temperature,  and  Hp=2kTp/mjVjn  where  Vjn  is  the  ion- 

neutral  collision  frequency.  For  vertical  field  lines 
sin^  1=1,  and  with  the  flux  S=NW£ji  (2.1)  becomes 


3N  _  3_ 

at  ~dz 


-NDa 

V 


n  a_N 
dz 


+ 


i  a Tn 
tp  az 


(q-l) 


(2.14) 


Expanding  (2.14)  we  get 


+  Da 


1 

9Tp 

dDa 

1 

DaTp 

9z 

dz 

'  Tp2 

[dzj 

1  92T  i 


9z2 


1  3Da  1  9Hq1 

*s*+ir-  ^2irJN  +  [^1] 


ip< 

(2.15) 


This  is  a  second  order,  parabolic,  partial  differential  equation  in  N(z,t)  of  the  form, 


-  A(z)  ^  +  B(z)  |j  +  C(z)  N  +  D(z) 
(2.16) 


The  coefficients  A,  B,  C,  and  D  are  altitude  dependent  parameters  which  can  be 
identified  with  the  quantities  in  brackets  in  (2.15). 

Analytic  solutions  are  possible  only  under  special  conditions,  e.g.,  diffusive 
equilibrium  (dN/dt=0)  and  simple  dependencies  for  the  coefficients.  Numerical 
solutions  permit  more  realistic  coefficients  and  even  allow  time  dependent 
coefficients  as  long  as  the  integration  time  step  is  sufficiently  short  for  them  to  be 
considered  constant.  This  became  crucial  after  we  added  a  calculation  for  Tc  which 
affects  the  plasma  temperature  Tp. 


A  diffusion  equation  requires  both  an  initial  condition,  in  this  case  an  altitude 
profile  of  0+  density  at  time  t=0,  and  boundary  conditions.  The  low  altitude  boundary 


is  chosen  low  enough  so  that  the  0+  density  calculated  from  local  chemistry  (CHEM) 
can  supply  one  condition.  At  the  upper  boundary  it  is  more  realistic  to  set  the  slope 
of  the  0+  density  which  is  related  to  the  flux  of  ions  across  the  boundary.  This  flux 
can  either  be  set  manually  or  estimated  from  scale  height  arguments.  Lacking 
measurements  or  reliable  estimates  of  the  flux,  one  can  test  for  a  range  of  values  into 
and  out  of  the  ionosphere.  There  is  a  critical  value  of  outgoing  flux  that  drives 
densities  negative  in  the  solution,  thereby  setting  an  upper  limit  of  upward  flux.  A 
similar  situation  also  occurs  in  the  Te  calculation. 

The  diffusion  equation  is  solved  numerically  by  a  Crank-Nicolson  scheme  (Gerald 
and  Wheatley.  1984)  starting  with  a  finite  difference  representation  of  (2.16). 


au/+1/2  4+1-u/ 


Al 

2 


i  +1 


at 

-  2u/  +  u 


i-  1 


i J+1 

i  + 1 


i+i  J+1 


-  2u r  +  u 


i- 1 


H2 


H2 


Bi 


J  -  J  J+1  -  J+1 

Ji  +1  ui  -  1  ui  +  1  ui  -  1 

—  +  — 


2H 


Ci  1"  j  j+1]  ui 

~K  +  ui  J  + 


2H 

d!  +  d/+1 


i 


(2.17) 


where  u(t,z)=u  is  density.  Index  j  corresponds  to  time  incremented  by  T  and  i 
corresponds  to  altitude  incremented  by  H.  Both  the  time  and  spatial  derivatives  are 
represented  by  central  difference  approximations  with  errors  of  order  Ax, 


1  3 


_  fi  +  1  ~  fi  - 1 
2Ax 

fi  +  1  -2fi  +  f 


2Ax 


+  O(Ax)2 
—  +  O(Ax)2 


Since  3u/3t  is  evaluated  at  the  midpoint  of  the  time  interval  (jJ+1)>  spatial 
derivatives  at  the  midpoint  are  approximated  by  averages  evaluated  at  the  beginning 
and  end  of  the  time  interval. 

Equation  (2.16)  expresses  unknown  densities  at  time  j+1  in  terms  of  known 
densities  at  t  .  Collecting  terms,  we  get 


2H2 
A: 


v 

4H 

£i_ 

PH2  '  4H 


H 


j+1 


( A:  C  ; 


_ l_ 

H2 

A 


H2  + 
(2.18) 


a  BJ) 

2  +  T  j 

Ci  B,' 
2  +  T 


♦ 


H+l 


_Al  V 

2H2  '  4H 


A: 


2H2  4H 


The  coefficients  in  (2.18)  are  determined  in  0+  DIFFUS  and  recast  in  the  following 
equation. 


Uji+f1  DMj  +  iji+1  EMj  +  FMj 
-  L|j_  1  DBj  +  iji  EBj  +  t|J+1FBj  +  Dj  =  BMj 

(2.19) 
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Within  the  main  loop  of  the  CHEM  program  the  coefficients  A,  B,  C,  D,  and  hence  DM  , 
EM  ,  etc.  are  constant  in  time.  The  right  side  of  (2.19)  contains  densities  at  the  end  of 
the  prior  time  step  j  (initial  altitude  profile  at  j=0,  t=0).  Since  (2.19)  couples  densities 
at  3  adjacent  altitudes  at  the  end  of  the  next  time  step,  j+1,  in  terms  of  quantities  from 
the  previous  time  step,  a  full  altitude  profile  at  time  tj  +  i  requires  solving  a  set  of 
linear  equations  for  u"j+!  i=0  to  n,  corresponding  to  z0  (low  altitude  boundary)  to  zn 

(high  altitude  boundary).  This  tridiagonal  system  is  solved  by  an  efficient  Gauss 
elimination  method  with  the  sole  restriction  that  the  altitude  spacing  H  be  less  than 
Aj/2Bj  for  all  i.  An  implicit  scheme  involving  several  altitudes  at  tj  +  j,  as  in  (2.19), 
has  the  virtue  of  stability  for  all  values  of  the  coefficients  as  opposed  to  explicit 
methods  which  handle  each  altitude  independently. 

Equation  (2.19)  applies  as  shown  except  at  the  boundaries,  zQ  and  zn  .  Since  the 
0+  density  at  the  lower  boundary  is  held  at  the  current  local  chemistry  value,  u^+1  = 

uJ  =  u  .  Thus  (2.19)  becomes 
o  o 

u^+1DMi  +  u!,  +  1EMi  +  uj>+1  FMi 
=  DBi  +  u!|  EBi  +  u;?  FBi 


or,  because  DM  =  -DB  , 


u^+1DMi  +uJ1  +  1EMi  +  uj>+1  FMi 
=  DBi  +  ul|  EBi  +  uj,  FBi 

(2.20) 
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To  incorporate  the  upper  boundary  flux  condition  into  the  solution,  we  begin  with 
(2.13)  (N=>u).  The  flux  of  0+  crossing  the  upper  boundary  (positive  upward)  is  given 
by 


4>n  =  unWd  = 


.  2  (1  dun  1  djr^ 

-unDa  s.n2|  I-  _  +  -  ^  + 


f  dun 


(2.21) 


A  positive  value  for  the  parameter  Qn  indicates  a  downward  net  flux  through  the 
upper  altitude  boundary.  From  (2.21),  the  slope  of  the  0+  density  at  the  upper 
boundary  can  be  specified  by  a  central  difference  approximation  that  gives  a  value 
for  un  +  l  required  by  (2.17). 


dun  Qn  -  Gun  un  +1  '  u  n  -  1 
dz  “  F  2H 


u 


+i  =  4i  - 1 


+  2H 


Qn  -  Gun 


Again  collecting  terms  as  in  (2.18), 


i  +  1  f  An^j  j +1  fAn  Cn'|  G  fAn  Bn)  1 
un-1  -  ^2  +  4  S2-  T  +  F  (Tr+  T  +  T 


j  fAn)  j  I”  1  (An  CnA  G  (An  Bn 

”  *"■  1  (h2 J  +  ^  [  T  ’  [&-  Tj-  -I_i  ~ 


F  H 


2fAn  Bn  'l 
CQ"  =  F  IT-  f  Qn 


+  CQn  +  Dn 


(2.22) 


At  the  upper  boundary  (2.19)  becomes 


Ujft1  DMn  +  ^+1EMn  + 

^  DBn  +  i||  EBn  +  CQn+  Dn  =  BMn 


(2.23) 


The  system  of  equations  represented  by  (2.19),  (2.20)  and  (2.21)  can  be  summarized  by 
a  n  x  n  matrix  equation. 


EMi  FMi  0 
DM2  EM2  FM2 

0  . 

0  DMj 
0  . 


EMi  FMi 

0  DMn-i 
0 


EMn-iFMn-i 
DMn  EMp 


(2.24) 


The  internal  time  step  dt=T  is  estimated  from  a  truncated  version  of  (2.14), 


To  evaluate  (2.25)  an  explicit  expression  for  Da  and  N  are  needed. 


Da 


Hpv 

g 


where  Hp  is  the  ion  (0+)  plasma  scale  height, 


H 


k(Te  +  Tj) 
mjg 


and  v  is  the  ion-neutral  collision  frequency.  The  e  and  i  indices  for 
and  mass  refer  to  electron  and  ion;  k  is  the  Boltzmann  coefficient  and 
gravitational  acceleration.  N(z)  can  be  approximated  in  terms  of  the 

„  Z~Z° 

N(z,t)  =N0(t)  e  HP 


leading  to, 

3N  _N_ 
dz  *  *  Hp 


the  temperatures 
g  is  the 
scale  height  by 


Substituting  these  into  (2.25)  we  get 


1 


H  p  d  D  a 
Da  dz 


=  Da 


N  f 

hp2  . 


o  r 


At 


AN  Hpg  ( 
N  n 


H  p  dDaY 1 
Da  dz  j 


(2.26) 


If  we  require  the  relative  density  change  to  be  small,  AN/N  =  e=(=.01,  e.g),  and  set 
Hpg/v  =x  ,  (2.26)  becomes 


_ et _  et 

\  H  p  dPa^|  ~  ( 

‘  Da  dz  J  ^  '  HdJ 

(2.27) 


The  second  term  in  the  denominator  of  (2.27),  the  ratio  of  the  plasma  scale  height  to 
the  diffusion  scale  height,  influences  the  final  5t.  At  high  altitudes  (600  km,  e.g.)  Hp 
is  typically  5-10  times  larger  than  Hd-  At  lower  altitudes  Hp/V  decreases  and  Hp=Hd, 
which  means  that  At  increases  rapidly.  Thus,  usually  only  the  highest  altitude  is 
required  to  determine  the  shortest  At  necessary  for  the  diffusion  calculation. 

The  following  examples  (Figures  2  and  3)  illustrate  some  diffusion  results.  They  all 
are  derived  from  the  same  incident  electron  flux  measured  on  the  PIIE  rocket  flight 
(Weber  et  al..  1987).  The  upper  and  lower  boundaries  arc  at  600  km  and  120  km, 
respectively;  the  altitude  spacing  is  10  km.  In  these  examples  both  the  electron  and 
ion  temperature  profiles  do  not  change  (Te=1000K,  Te=2000K,  TC=1500K,  at  600  km). 

The  effect  of  diffusion  can  be  seen  most  readily  by  comparing  Figure  2  (local 
chemistry  only)  with  Figure  3  (local  chemistry  and  0+  diffusion).  Both  figures  begin 


with  the  same  initial  electron  (=total  ion)  density  profile  at  t=0,  and  show  successive 
electron  density  profiles  at  10  min  intervals  up  to  1  hr.  With  local  chemistry  acting 
alone,  rapid  production  at  300  km  from  the  relatively  soft  incident  electron  flux 
causes  a  rapid  buildup  of  ionization  above  250  km.  Below  200  km  equilibrium 
densities  are  reached  after  10  min.  Diffusion  essentially  drives  0+  ions  away  from 
steep  density  gradients,  flattening  the  sharp  density  peak  produced  by  local 
chemistry.  As  expected  the  chemistry  and  diffusion  results  are  virtually  identical  at 
low  altitudes  (<200  km).  Higher  up,  the  density  approaches  equilibrium  more  rapidly 
with  diffusion;  the  peak  density  is  also  slightly  lower  in  altitude  and  only  about  50% 
of  the  chemistry  result.  However,  because  of  the  flatter  gradient  ,  the  O+  density 
after  1  hr  at  600  km  is  3.2  x  104  (cm'3)  compared  to  3.5  x  10^  without  diffusion. 

Figures  4a-d  show  the  effect  of  varying  the  high  altitude  0+  flux  from  1  x  10^  m*2 
s'1  (down)  to  -5  x  lO1^  (up)  The  differences  are  hardly  noticeable  undl  the  critical 

outflow  value  is  neared  (Figure  4d),  which  has  decreasing  density  above  400  km. 

After  10  min  the  diffusion  speed  is  100  m/s  down  in  Figure  4a  and  2  km/s  up  in 
Figure  4d. 

D.  Ionospheric  Temperature 

The  thermal  energy  balance  for  ionospheric  electrons  is  represented  by  a  one 
dimensional,  diffusion  type  equation  (Hastings  and  Roble.  1977) 


aE  3  ,3Te  1  B  (  „  BJV) 

a7=  2  nek"at  =a  all3  Ke  if) 


+  Q-L 


(2.28) 


E  =  energy  density  of  thermal  electrons  (eV  cm" 3) 
nc  =  electron  density  ( cm" 3) 
k  =  Boltzmann  constant 


20 


Ke  =  thermal  conductivity  of  the  electron  plasma  (eV  ernes'1  K'1) 

a  =  cross  sectional  area  of  magnetic  flux  tube 

Q,  L  =  energy  gain  and  loss  rates,  respectively  (eV  cm‘3  s'1) 

The  same  Crank-Nicolson  finite  differencing  scheme  that  solves  for  0+  diffusion  is 

used  to  solve  the  heat  conduction  equations.  However,  0+  diffusion  is  relatively 
straightforward  since  the  coefficients  depend  only  on  temperature  and  other 

parameters  that  are  updated  at  regular  time  interval.  Temperature,  on  the  other 

hand,  requires  solving  a  non-linear  equation  with  coefficients  that  depend  on  the 

temperature.  With  short  time  steps  the  problem  is  managed  by  assuring  that 
coefficients  change  slowly  enough  to  be  linearly  projected  to  the  end  of  each  time 
step.  Alternatively,  Hastings  and  Roble  (1977)  and  Mantas  et  al.  (1981)  use  Newton's 
method  to  solve  a  different  formulation  of  the  non-linear  diffusion  type,  heat 
conduction  equation. 

Expanding  (2.28) 


d  a 
dz 


K< 


L 


(2.29) 


b  =  2/3nek. 

In  the  form  of  (2.16)  the  coefficients  are, 

dzj  ’  C  =  0 .  D  =  Q  -  L  . 

The  thermal  conductivity  given  by  Banks  and  Kockarts  (1973),  Vallance-Jones  (1974) 
and  others  is  expressed  as 


A  =  bKe  ,  B  =  b 


K  7.7x1  0s  Te5/2 

1  +  3. 22x10*  ^XnmQDk 
n©  m 


2  I 


Qdr-4  /fkiyvwl  (ve(Ng>  +  Ve(°2)  +  Ve(0)) 

3  I^Tcme  > 

nk  =  density  of  neutral  species  N2,  O2,  0  (cm‘3). 

QDk  involves  collision  frequencies  of  electrons  with  the  neutral  constituents,  which 
in  turn  depend  on  Te  and  nfc.  The  derivative  of  the  flux  tube  area  is  calculated  by 
assuming  r’3  dependence  for  the  cross  sectional  area  of  a  dipole  flux  tube.  In 
practice  this  modification  turns  out  to  be  quite  small,  but  it  is  retained  for 
completeness. 

The  final  term  in  (2.29)  represents  the  net  external  source  of  heat  given  by  the 
difference  between  heat  generated  by  the  electron  flux  and  losses  to  the  generally 
colder  neutrals  and  ions.  A  normalized  heat  input  from  the  electron  flux  transport 
code  is  multiplied  by  ne  to  yield  Q.  Losses  tabulated  in  Vallance-Jones  (1974)  are 
categorized  as  cooling  by  elastic  collisions  with  neutrals  and  ions,  cooling  by 
rotational  and  vibrational  excitations  of  N2  and  O2,  and  cooling  by  electronic 
excitations  of  O. 

Within  the  program  unit  THERM  in  Figure  1  ,  which  shows  the  overall  program 
logic,  there  are  several  subroutines  needed  to  perform  the  temperature  calculation. 
SETKE,  which  is  executed  only  at  the  beginning  of  the  program,  sets  up  the  altitude 
grid  (now  in  10  km  steps  from  130  to  600  km  altitude)  and  initializes  the  temperature 
profile.  If  one  selects  the  option  to  calculate  electron  temperature,  the  controlling 
routine  ETEMP  is  called  immediately  after  some  input  and  initialization  routines  are 
run . 

Because  Te  can  change  rapidly  from  its  initial  profile,  depending  on  conditions, 
the  first  times  steps  should  be  short  (.05  -  .1  s,  e.g.).  Temperature  and  densities  at  the 
beginning  of  the  time  step  are  used  to  compute  coefficients  in  TECOEFF,  which  is 
called  from  the  numerical  solution  routine,  CONDTE.  In  a  major  departure  from  the 
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diffusion  routine,  coefficients  values  at  the  end  of  the  time  step  are  linearly 

extrapolated  from  the  two  previous  time  steps  (except  at  t=0)  to  reflect  potentially 

large  changes  in  the  coefficients.  A  quadratic  extrapolation  was  also  tried  but  the 

results  did  not  differ  appreciably  from  the  linear  version.  Again,  at  the  upper 

boundary  the  user  has  a  choice  between  specifying  a  temperature  or  heat  flux  (Qn  = 
9Tn 

-Kn“ —  ,  n=upper  altitude  index). 

a  z 

At  the  end  of  a  time  step  DTTH,  CONDTE  returns  control  to  ETEMP,  which,  after  the 
first  CONDTE  call,  reruns  CONDTE  with  new  time  steps,  dt=DTTH/2  and  dt=2  x  DTTH. 
Results  of  this  testing  procedure  are  printed  at  the  selected  output  times  and  allow  the 
user  to  decide  whether  DTTH  should  be  shortened  or  lengthened.  Normally,  ETEMP 
executes  CONDTE  until  the  end  of  an  update  interval,  DTUP,  is  reached.  Typical  values 
for  DTUP  used  in  the  examples  range  from  .Is  for  t<ls  to  5s  for  t>30s.  If  the  ion 
temperature  option  is  also  selected,  similar  routines  are  executed  to  calculate  the  ion 
temperature. 

After  the  temperatures  are  calculated,  0+  diffusion  coefficients  and  chemical  rate 
coefficients  are  recalculated  using  the  new  temperature  profiles.  The  local 
chemistry  and  diffusion  routines  then  run  to  find  new  densities  at  the  end  of  the 
DTUP  interval,  after  which  the  main  program  checks  whether  it  is  time  to  output 
results  or  end  the  model  run.  If  neither  condition  is  satisfied  the  program  loops  back 
for  a  new  temperature  calculation  to  the  end  of  the  next  DTUP  interval  and  so  on. 
Because  temperatures  are  computed  first,  rate  coefficients  for  chemistry  (constant 
throughout  the  DTUP  interval)  are  based  on  temperatures  at  the  end  of  the  interval; 
thus  the  temperature  is  "ahead"  of  the  density  by  the  DTUP  interval.  Some  form  of 
averaging  could  minimize  the  discrepancy,  but  by  keeping  the  DTUP  interval 
relatively  short  compared  to  the  time  scale  of  density  changes,  .1-1  Os  compared  to 
several  minutes  in  most  cases,  the  mismatch  should  be  negligible. 
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Because  of  the  complexity  of  the  model  it  is  very  difficult  to  make  direct 
comparisons  with  measured  temperatures.  Although  most  results  appear  plausible,  a 
real  situation  is  usually  much  more  dynamic  than  the  constant  flux  assumed  in  the 
model;  moreover,  temperatures  respond  much  more  quickly  to  changes  in  flux  input 
than  densities,  particularly  in  the  F-region.  Like  diffusion,  the  high  altitude  heat 
flux  is  generally  not  available,  although  a  high  altitude  temperature  can  be  specified 
from  measurements  if  desired. 

Even  though  an  analytic  solution  to  the  heat  conduction  equation  with  variable 
coefficients  is  not  possible,  a  steady  state  solution  with  constant  (in  time)  coefficients 
can  be  expressed  as. 


T'z>  -  Ti  -  tf) d2' +  £  A^hS(z")dz'>2' 

(2.30) 

z  =  altitude,  zi  =  low  altitude  boundary,  zh  =  high  altitude  boundary 

3T 

S(z)  =  -(Q(z)-L(z)],  ph  =  T~  at  zh. 

dz 


A  limited  check  of  the  Crank-Nicolson  solution  used  in  the  model  was  carried  out  by 
stripping  the  temperature  routines  from  the  full  code,  running  the  reduced  model  to 
steady  state,  and  comparing  the  results  with  a  numerical  integration  (Simpson  rule) 
of  (2.30).  Table  2  shows  the  results  of  one  comparison.  The  table  indicates  that  the 
model  starts  with  Te  at  about  2050  K  at  600  km  altitude,  which  then  cools  to  about  1500 
K  after  960  s  (ionosphere  time).  Analytic  results  using  the  same  coefficients  are 
nearly  identical  to  the  model  results,  giving  some  confidence  that  at  least  the 
numerical  procedure  is  correct  for  a  special  situation  with  realistic  temperatures. 
Although  this  check  is  a  necessary  condition  for  the  validity  of  the  model,  it 
certainly  is  not  a  conclusive  test  under  more  realistic  conditions. 
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3.  Running  the  Ionospheric  Chemistry  Model 

The  growing  complexity  of  the  chemistry  codes  prompted  the  creation  of 
setup  programs  that  allow  a  knowledgeable  user  to  more  easily  specify  the 
numerous  files  and  parameters  required  by  the  model  program  itself.  Even 
with  this  aid  the  model  cannot  be  called  "user  friendly",  but  the  short 
preparatory  codes  and  procedures  do  reduce  the  likelihood  of  typing  errors 
and  missing  input  files,  permitting  faster  turn  around  especially  when 
submitting  similar  jobs  with  varying  input  parameters.  The  various  options 
and  menus  available  are  described  in  detail  below. 

The  setup  program  runs  in  two  parts;  part  1  asks  for  input  file  names, 
output  file  names,  and  some  processing  options;  part  2  allows  one  to  change  or 
accept  parameters  controlling  boundary  conditions,  initial  conditions,  and 
output  options.  Once  the  selections  are  made,  the  ionospheric  model  can  be 
either  submitted  as  a  batch  job  or  run  interactively  at  the  terminal. 

To  run  the  programs  a  reasonably  good  working  knowledge  of  the  Cyber 
NOS-2  operating  system  is  desirable  and  almost  essential,  particularly  when 
problems  occur,  which  inevitably  happens.  Procedures  are  used  to  automate 
the  model  as  much  as  possible;  these  also  require  several  indirect  and  direct 
access  permanent  files.  The  modified  model  codes  arc  still  written  in  FORTRAN 
version  4,  as  were  the  original  codes.  Since  this  version  will  eventually  be 
unsupported  they  should  be  converted  to  the  new  standard,  version  5.  A 
conversion  has  been  done  but  the  compiled  code  took  substantially  more  CPU 
time  for  execution  in  some  cases.  This  problem  may  be  connected  with 
printing  options,  DO  loop  procedures  or  even  the  Post  Mortem  Dump  utility, 
and  until  the  exact  solution  is  found  it  was  decided  to  stay  with  a  relatively 
efficient  working  code.  The  original  documentation  was  adequate  but  spotty; 
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additional  comments  have  been  added  at  crucial  points  in  the  old  code  as  well 
in  the  new  or  modified  portions. 

The  chemistry  model  is  started  by  typing  and  entering  (RETURN)  the  Cyber 
procedure  CHEM  (on  a  permanent  file  with  same  name).  A  list  of  options  then 
appears  with  brief  descriptions.  When  satisfied  with  the  option  list,  type  and 
enter  GO  to  begin  a  second  procedure  (BCHEMCZ)  which  contains  the  editing 
program  for  all  model  parameters.  One  can  choose  to  edit  a  single  group  of 
parameters  or  the  entire  list.  Depending  on  the  selection,  the  current  value  of 
a  parameter  is  displayed  with  the  option  to  change  or  not;  it  is  possible  to  go 
back  and  reedit  any  selection.  Commonly  used  parameter  sets  can  be  saved  on 
an  indirect  access  permanent  file  which  is  specified  as  an  option  on  the  first 
setup  program  (FORMATTED  FILE  FOR  CHEM-).  The  last  used  parameter  set  is 
also  saved  on  a  local  file  named  PREV.  Again,  when  satisfied  with  the 
parameter  settings,  type  GO  to  submit  the  model  as  a  batch  or  interactive  job. 

A.  First  List  (Options,  LIST=PART) 

After  typing  and  entering  CHEM  a  list  of  options  with  the  current  settings 
appears;  the  choices  can  be  grouped  in  the  categories  of  input  file,  output 
files,  and  processing  options. 

1.  Input  files 

PRFILE  -  Production  rate  file  created  by  the  electron  flux  transport  (B3C) 
and  production  rate  codes  (PRATES)  executed  in  the  TRAN  procedure.  This  file 
contains  the  altitude  grid,  neutral  atmosphere  model,  production  rate  profiles, 
and  the  normalized  electron  plasma  heating  rate;  enter  this  file  name  without 
the  PR  prefix. 

INITIAL  ION  DENSITY  -  A  deck  name  in  the  DATAPL  (UPDATE  utility)  file 
denoting  a  file  that  contains  specified  initial  density  profiles  of  selected  ions. 
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NEUTRAL  ATM.  MODEL  -  A  deck  name  in  the  DATAPL  file  that  contains  the 


appropriate  neutral  atmosphere  model;  this  is  not  used  if  the  neutral 
atmosphere  model  from  the  PRFILE  is  specified  from  the  parameter  list 
(below). 

ALTERNATE  INITIAL  DEN  FILE  -  A  file  (OUTPUT  FILE  (ALL  DENSITIES)) 
created  by  a  previous  chemistry  run,  containing  density  and  temperature 
profiles  at  the  last  output  time,  that  can  be  used  in  place  of  the  file  designated 
by  INITIAL  ION  DENSITY.  This  is  useful  if  initial  densities  based  on  a  particular 
flux  are  desired,  for  example  a  relatively  steady  flux  that  produces  a 
background  ionosphere. 

2.  Output  files 

PLOTFILE  (DEN.,  TEMP.,  ELEC.  COND.)  -  A  file,  used  by  plotting  programs,  that 
contains  calculated  electron  and  0+  density  profiles,  as  well  as  electron 
temperature,  ion  temperature,  and  the  Hall  and  Pedersen  electric 
conductivities.  Data  is  collected  at  output  times  specified  by  the  user  in  the 
parameter  list.  A  FORTRAN  plotting  routine  on  the  Cyber  has  been  developed 
(PLTESTP  for  pen  plots),  or,  for  interactive  plots  on  the  screen,  there  are 
ZB  ASIC  programs  (with  special  IGL  plotting  commands)  for  the  Zenith  Z-100 
personal  computer. 

OUTPUT  FILE  (ALL  DENSITIES)  -  A  file  containing  all  density  and 
temperature  profiles  generated  at  the  last  output  time  for  the  run.  This  is  the 
ALTERNATE  INITIAL  DEN  RLE  specified  by  a  later  run  if  one  wished  to  use 
these  profiles  as  initial  conditions. 

?  F  occssing  options 

(B)ATCH  OR  (I)INTERACTIVE  JOB  -  Enter  B  for  a  batch  job  that  is  routed  to  the 
input  queue  for  execution;  enter  I  for  immediate  interactive  execution. 
Interactive  jobs  are  useful  for  testing  because  execution  usually  starts 
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immediately,  although  the  users  terminal  session  is  tied  up  until  completion 
(unless  the  job  is  detached).  Batch  processing  is  generally  desirable  for 
routine  runs  because  multiple  runs  can  be  submitted  and  the  terminal  is  freed 
up. 

PARAMETER  FILE  -  a  file  containing  the  parameters  that  can  be  edited 
(below).  Normally  this  is  the  DCPL  file  associated  with  the  BCHEMED 
procedure;  however,  another  file  in  the  same  format  can  be  substituted. 

BATCH  ROUTE  CODE  (IN, TO)  -  a  parameter  designating  whether  the  printed 
model  output  appears  on  the  OUTPUT  (select  IN)  file  printed  at  the  computer 
center  or  on  a  wait  file  which  can  be  accessed  at  the  terminal  by  a  QGET 
command  (select  TO).  TO  is  the  normal  choice  since  the  user  can  quickly  check 
the  results,  which  can  be  ROUTED  to  a  line  printer  using  the  DC=PR  option.  The 
output  file  is  usually  to  long  to  be  printed  locally,  but  it  can  be  edited  so  that 
portions  may  be  printed,  if  desired. 

SEQUENCE  NO.  (DEFAULT=1)  -  If  multiple  batch  jobs  are  submitted  each  run 
must  have  a  distinct  sequence  number  so  that  each  will  access  the  correct 
files. 

The  following  items  are  listed  only  if  the  LIST=ALL  parameter  is  explicitly 
set  when  the  CHEM  procedure  is  invoked.  Most  are  relevant  only  for 
development  purposes  and  should  be  changed  only  if  the  user  is  familiar  with 
the  programs  and  procedures. 

RATE  COEFF  CORRECTION  FILE  -  an  UPDATE  correction  file  that  contains 
modified  and  additional  chemical  reactions  pertinent  to  ionospheric  chemistry 
(collected  from  Vallance-Jones.  1977).  The  default  setting  is  PICORR. 

FTH  (THERMAL  CODE  SUFFIX)  -  a  character  used  to  specify  the  file  containing 
the  temperature  subroutines.  The  default  setting  is  4. 
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FEL  (ELEC.  COND.  CODE  SUFFIX)  -  a  character  used  to  specify  the  File 
containing  the  electric  conductivity  subroutines.  The  default  setting  is  Y. 

PRFILE  FROM  OLD(O)  OR  NEW(l)  CODE  -  choice  determines  whether  PRFILE 
format  is  that  generated  by  original  FORTRAN  4  code  or  revised  FORTRAN  5 
version.  The  default  is  0,  designating  the  FORTRAN  4  version. 

NEW  COMPILATION  (Y,N)  -  choice  determines  whether  previously  compiled 
model  code  is  executed  (N),  or  new  source  code  is  to  be  compiled  before 
execution.  The  default  setting  is  N. 

OPTICS  (Y,N)  -  choice  determines  whether  File  containing  source  code  for 
calculation  of  optical  emissions  is  merged  with  main  chemistry  code  before 
compilation  (relevant  only  if  Y  is  selected  for  NEW  COMPILATION  option, 
above).  The  default  setting  is  Y. 

CHEM  PROCEDURE  VERSION  -  a  character  sufFix  specifying  the  File 
containing  the  procedure  that  executes  the  ionospheric  chemistry  model.  The 
default  setting  is  4. 

PROGRAM  STOP  PARAMETER  -  a  number  allowing  a  programmer  to  stop 
model  execution  at  a  desired  point  in  the  code  by  inserting  STOP  statements  in  a 
modiFied  sor.rce  code.  The  default  setting  is  0. 

B.  Second  list 

After  the  editing  program  is  loaded  the  main  list  appears: 

A.  INITIAL  DENSITY  FILES 

B.  MAIN  CHEM  OPTIONS 

C.  DIFFUSION  ALT.  GRID 

D.  TE,TI  UPPER  ALT.  FLUXES 

E.  CHEM  TIME  STEPS 

F.  SCALE  FACTORS 

G.  0+  UPPER  ALT.  FLUX 
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H.  DIFFUSION, TEMP.  PARAMETERS  AND  OUTPUT  OPTIONS 

I.  DISPLAY  ALL  CHOICES  (A-I) 

ENTER  SELECTION  A-I 

OR  TYPE  GO  TO  RUN  WITH  CURRENT  PARAMETERS 
After  making  a  selection  a  group  of  parameters  (from  the  DCPL  file  specified 
in  the  first  list)  is  displayed  along  with  the  option  of  making  a  change  or 
letting  them  stand.  When  satisfied  with  the  displayed  parameters  the  user  can 
return  to  the  main  list  and  choose  another  group  to  edit  or  enter  GO  to  submit 
the  job  for  execution.  Usually  after  an  incorrect  response  the  CHANGE(Y/N)? 
prompt  reappears.  If  a  mistake  is  made  in  entering  data  it  can  be  edited  again, 
or  if  hopelessly  lost,  the  CONTROL-T  abort  combination  can  be  entered  to  leave 
the  procedure  and  return  to  the  NOS  backslash  prompt  (\),  after  which  the 
CHEM  procedure  can  be  tried  again.  After  entering  GO  from  the  main  list  or 
the  H  option,  the  status  report  INPUT  FILE  WRITTEN  appears  along  with  a 
listing  of  the  batch  job  file  (if  selected)  that  will  be  routed.  The  job  is  then 
routed  and  the  status  after  submission  (INPUT  QUEUE  or  EXECUTING,  e  g.)  is 
listed  on  the  screen.  If  the  interactive  option  is  chosen  from  the  CHEM 
argument  list  the  job  will  run  until  the  EXIT  response  signals  the  job  has 
finished.  While  the  job  is  running  one  can  enter  ESC(APE)  and  then  enter  E  to 
get  a  status  report.  Model  output  from  an  interactive  job  will  be  on  the  local 
file  OUTPR. 
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4.  Parameterization  of  the  Flux  Transport  Model 


A.  Motivation 

Besides  the  inherent  complexity  of  the  ionospheric  chemistry  model,  the 
code  requires  many  parameters  and  external  inputs.  For  example,  there  arc 
initial  and  boundary  conditions,  a  background  neutral  atmosphere,  and 
production  of  various  atomic  and  molecular  species  as  well  heat  generated  by 
incident  electrons,  protons,  or  solar  radiation.  In  many  situations  these 
parameters  do  not  change,  but  particle  fluxes  are  often  highly  dynamic  and 
spatially  nonuniform,  especially  near  auroral  arcs.  Since  the  models 
(production  rates  and  chemistry)  are  one  dimensional  along  a  magnetic  field 
line,  they  do  not  account  for  horizontal  transport  produced  by  electric  fields. 
However,  the  models  could  be  run  with  fluxes  measured  at  various  points 
across  the  arc,  and,  by  combining  the  measured  electric  field  with  the 
calculated  horizontal  density  gradient,  an  additional  source  term  could  be 
included  in  the  chemical  rate  equations.  Similarly,  a  time  dependent  flux 
could  be  simulated  by  inserting  different  production  rates  at  appropriate  times 
in  the  chemistry  calculations. 

The  calculation  of  production  rate  profiles  by  an  incident  particle  flux 
alone  is  a  complicated  task.  There  are  semi-empirical,  range  theoretic  models 
(e  g-,  Rees  (1963;  Roble  and  Rees  (1977))  and  quantitative  models  based  on  the 
Boltzmann  equation  and  cross  sections  for  ionization  and  excitation 
(Strickland.  1976;  Banks  et  al..  1974;  Mantas  et  al..  1975).  We  have  applied  the 
Strickland  model  to  studies  of  polar  cap  auroral  arcs  (Weber  et  al..  1988). 
Production  rates  calculated  for  the  maximum  electron  flux  measured  by  a 
rocket  over  an  arc  were  input  to  the  ionospheric  chemistry  model  described 
above.  Since  the  models  (flux  transport  and  chemistry)  are  one  dimensional 
along  a  magnetic  field  line,  they  do  not  account  for  horizontal  transport 


produced  by  electric  fields.  However,  the  models  could  be  run  with  fluxes 
measured  at  various  points  across  the  arc,  and,  by  combining  the  measured 
electric  field  with  the  calculated  horizontal  density  gradient,  an  additional 
source  term  caused  by  E  x  B  convection  could  be  included  in  the  chemical  rate 
equations. 

A  need  for  production  rate  profiles  corresponding  to  many  different  fluxes 
is  restrictive  because  of  the  large  computational  expense  involved.  Aside  from 
multiplication  by  a  constant  factor  or  a  simple  time  dependence,  there  has 
been  no  alternative  to  a  new  run  of  the  flux  transport  code  for  each  unique 
flux.  We  have  found,  however,  that  production  rate  profiles  can  be  quickly 
approximated  by  scaling  a  table  of  rate  profiles  calculated  for  a  specified  unit 
flux  over  a  grid  of  discrete  energies.  We  will  describe  below  how  the  tables  are 
generated  along  with  several  examples. 

B.  Outline  of  Strickland  Flux  Transport  Model 
The  Boltzmann  equation  for  electrons 


a t 

at 
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(4.1) 


is  reduced  to  an  equation  of  transfer  with  azimuthal  symmetry  about 
geomagnetic  field  lines, 


a<j>(t,E,fl) 
4  al 


=  -  <t)(x,E,}i)  +  I 
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0  =  flux  (cm'2  s*1  sr1  eV*1),  R]  =  redistribution  function 
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=  cosine  pitch  angle,  k  corresponds  to  neutral  constituent  species  (N2,  O2, 


O) 

X  =  scattering  depth,  dt  =  Znk(z)ok{E)dz,  n^=  neutral  density,  z  =  altitude 

=  total  cross  section  for  all  scattering  processes 
rk  =  nk(z)ok(E)/In.(z)oi(E) 


This  form  of  the  Boltzmann  equation  is  one  dimensional  and  neglects  any 
forces,  e.g.,  electric  fields,  along  magnetic  field  lines.  This  permits  the  model 
to  step  successively  from  higher  to  lower  energies  since  an  electron  can  only 
change  direction  by  elastic  scattering  or  lose  energy  by  ionization  or 
excitations  of  the  neutral  constituents. 

A  matrix  representation  of  the  equation  of  transfer  is  then  derived  over  a  2- 
dimcnsional  grid  of  energies  and  pitch  angles.  The  coefficients  arc  obtained 
by  assuming  that  the  flux  within  each  grid  cell  can  be  approximated  by  a 
function  quadratic  in  log(E  =  energy)  and  linear  in  .  A  power  law  continues 

the  flux  beyond  the  highest  grid  energy.  At  each  energy  the  flux  is  given  by 


90i(x,E,4) 
w  at 


+  Sj' 


i,j  indices  range  over  pitch  angle;  n,m  range  over  energy. 

Sj'  =  Er.  Z  A  <J>mj  =  source  function  from  higher  energies, 
k  K  jnem  nmij 

An  eigenvalue  solution  of  the  matrix  equation  is  calculated  at  an  intermediate 
altitude  of  the  altitude  range.  The  high  altitude  incident  flux  constitutes  one 
boundary  condition;  zero  upward  flux  at  the  lowest  altitude  completes  the 
necessary  houndary  conditions.  Differential  flux  solutions  at  other  altitudes 
are  approximated  by  an  iterative  method.  Finally,  the  computed  profile  of 


differential  electron  flux  is  integrated  over  pitch  angle  and  input  to  a 
production  rate  code  that  yields  production  rate  profiles  of  14  ionic  and 
neutral  species. 

The  principal  requirements  of  the  model  are: 

(1)  Tables  of  ionization  and  excitation  cross  sections  for  electrons  colliding 
with  N2,  O2,  and  O. 

(2)  A  neutral  atmosphere  model  with  densities  of  N2,  O2,  and  O.  The  examples 
below  use  a  Jacchia  1977  model  having  1000  K  exospheric  temperature. 

(3)  Specified  energy  and  pitch  angle  grids.  Since  we  are  mainly  concerned 
with  soft  electron  fluxes  often  observed  over  polar  cap  auroral  arcs,  typically 
having  energies  around  a  few  hundred  eV,  we  use  19  transport  energies 
ranging  from  1800  eV  to  10  eV  and  7  energies  from  8  eV  to  3  eV  which  are 
deposited  locally.  20  pitch  angles,  10  each  over  the  up-  and  downgoing 
hemispheres  are  standard. 

(4)  A  differential  flux  at  the  grid  energies  and  pitch  angles  incident  at  the 
upper  boundary.  The  examples  assume  a  downgoing  isotropic  flux  at  450  km 
altitude,  near  the  altitude  of  rocket  observations  in  recent  AFGL  experiments 
(PIIE,  Polar  Arcs). 

C.  Constructing  the  Tables 

When  results  for  different  incident  fluxes  incident  on  the  same  neutral 
atmosphere  are  desired,  computations  involving  (1),  (2),  and  (3)  above  are 
duplicated.  Because  each  transport  run  needed  to  calculate  a  flux  profile  can 
take  nearly  2  min  CPU  on  a  Cyber  760  computer,  situations  where  dozens  of 
runs  arc  required  may  consume  excessive  amounts  of  machine  time.  Execution 
time  can  be  shortened  considerably  by  determining  the  contribution  to 
production  by  a  unit  flux  at  each  grid  energy.  Multiplying  each  contribution 
by  the  actual  flux  at  that  energy  scaled  to  the  unit  flux  yields  the  contribution 
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by  that  flux.  Final  rate  profiles  are  produced  by  adding  the  contributions  from 
all  energies.  A  complete  set  of  rate  profiles  from  a  single  differential  flux 

spectrum  can  be  generated  in  1  s  of  CPU  time. 

Although  ultimately  there  may  be  better  ways  to  do  the  job,  Figure  5 
illustrates  how  to  construct  the  tables.  This  schematic  shows  how  the 
contributions  from  energy  E5  is  found.  Spectrum  4  (S4)  represents  a  transport 
and  production  run  yielding  production  rate  profiles  F4  for  a  constant 
incident  flux  <j)  1  up  to  E4  and  a  lower  flux  <j>o  at  higher  energies  on  the  grid.  <))  1 
is  comparable  to  the  level  of  a  realistic  flux,  <J>o  is  low  enough  to  add  essentially 
no  production  compared  to  <J)j.  In  the  examples  below  <J>]  =1  x  10^ ,  <}>o  =  *  x  10®.  S5 
is  similar  to  S4  but  <j>i  extends  to  E5.  By  subtracting  the  production  rate  profiles 
generated  by  S4  from  the  corresponding  profiles  generated  S5,  one  is  left  with 
the  production  by  <|)  1  at  E5.  Subtracting  neighboring  pairs  of  runs  at  the  other 

grid  energies  completer  the  table  of  contributions  from  all  energies  except 
one  (  the  lowest  transporting  energy). 

A  range  of  20  energies,  for  example,  still  requires  20  full  transport  runs  to 
construct  a  full  table.  This  also  demands  considerable  computer  time;  the 
tables  would  be  justified  only  when  a  relatively  large  number  of  fluxes  is 
encountered,  such  as  when  modelling  detailed  behavior  across  an  auroral  arc 
or  in  global  3-D  ionospheric  simulations.  Global  modelling  would  require  a 

table  for  each  neutral  atmosphere,  but  careful  interpolation  may  lessen  the 
need  for  many  such  table.  Similarly,  each  pitch  angle  dependence,  isotropic, 
beamlike,  or  conic  would  demand  its  own  table. 

The  time  to  generate  a  table  can  be  further  shortened  by  considering  some 
details  of  the  code  logic.  Namely,  (1)  eigenvalues  need  not  be  recalculated  for 
each  run,  and  (2)  source  functions  and  flux  results  at  higher  energies  can  be 


saved  and  reused  Tor  subsequent  runs  at  lower  energies.  These  savings  could 
reduce  the  time  to  construct  tables  by  50%. 

D.  Results  and  Examples 

Figures  6a-f  are  plots  of  tabular  rate  profiles  of  3  ionic  species  for  selected 
energies  ranging  from  55  to  1300  eV  in  response  to  isotropic  incident  flux  of 
electrons,  <j>i  =  1  x  105  (cm'2  s'1  sr1  eV'1)  at  each  energy.  The  altitude  of  peak 

production  ranges  from  over  300  km  at  55  eV  to  150  km  at  1300  eV;  peak  rates 
range  from  10' 1  to  102  (cm'3  s'1).  Note  that  the  production  rate  scale  changes 
at  higher  energies. 

To  compare  table  results  with  a  full  electron  transport  and  production  run, 
a  flux  spectrum  from  the  AFGL  PIIE  experiments.  Figure  7a,  was  input  to  the 
full  code.  Figure  7b  shows  production  rates  from  the  tables  using  the  flux  in 
Figure  7a.  Table  3  compares  rates  scaled  from  the  tables  with  those  from  the 
full  code  (shown  in  parentheses  when  different).  It  is  evident  that  most  table 
results  are  very  close  to  the  full  model  values.  Profiles  generated  from  the 
tables  for  a  200  eV  Maxwellian  flux  (Figure  8a)  are  plotted  in  Figure  4b.  The 
type  of  agreement  is  similar  to  that  seen  in  Table  3. 

Work  remains  to  explore  where  and  why  discrepancies  exist  and  to  place 
limits  on  the  types  of  fluxes  the  method  is  able  to  handle.  The  examples  above 
and  Maxwellian  fluxes  at  80  and  600  eV  (not  shown)  indicate  that  the  tables 
perform  well  for  reasonably  behaved  fluxes.  Further  tests  are  needed  for 
higher  energy  auroral  precipitation  at  several  keV  and  for  more  unusual 
fluxes,  e.g.,  double  peaks.  A  cursory  glance  at  Table  3  and  the  other  species 
suggests  that  the  ionic  production  rates  compare  better  than  some  neutral 
excited  species.  Thus,  the  table  may  be  better  suited  for  studies  involving  total 
ion  density  rather  than  excited  states. 
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As  indicated  earlier,  auroral  arc  modelling  could  be  an  immediate 
beneficiary  of  the  tabular  production  rates.  The  inherent  one  dimensional 
nature  of  the  transport  model,  however,  assumes  horizontal  uniformity  and 
limits  the  resolution  of  a  structure  to  scales  greater  than  an  electron 
gyroradius.  This  is  not  a  serious  concern  for  arcs  several  tens  of  kilometers 
across.  The  tables  will  permit  much  more  experimentation  with  spatial 

structure  and  time  dependence  to  simulate  the  increasingly  detailed 
observations  provided  by  radars  and  photometers. 

Finally,  a  longer  term  goal  would  be  to  set  up  a  data  base  of  lookup  tables  for 
users  to  access  production  rates  for  a  wide  variety  of  neutral  atmospheres, 
incident  altitudes,  and  pitch  angle  dependences.  Several  altitudes  may  be 
required  to  accommodate  rocket,  satellite,  or  magnetospheric  model  inputs. 
Global  modelling  demands  the  widest  range  of  possibilities,  but  could  gain 
considerably  by  refining  the  cruder  zonal  inputs  now  used. 
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5.  PIIE  Results 


Figure  9  displays  measurements  taken  by  the  instrumented  rocket  payload 
as  it  passed  over  a  complicated  system  of  F-layer  auroral  arcs  near  apogee.  The 
vertical  line  indicates  the  time  chosen  for  model  comparisons.  The  electron 
density  had  a  local  maximum  of  3  x  104  cm'3  and  Te  was  about  3000  K  at  415  km 
altitude.  Electrons  precipitating  into  the  arc  carried  an  energy  flux  of  about  .1 
(erg  cm'3  s'3),  nearly  an  order  of  magnitude  greater  than  in  adjacent  regions. 
A  Jacchia  1977  neutral  atmosphere  model  with  1000  K  exospheric  temperature 
was  employed.  Densities  at  450  km  and  105  km  altitude  for  N2,  O2,  and  O  range 
from  1.3  x  106  to  3.9  x  1012  cm’3,  2.9  x  104  to  7.7  x  1011  cm'3,  and  4.7  x  107  to  3.2 
x  103  3  cm-3,  respectively.  Energy  transfer  from  the  energetic  electrons  to 
ambient  thermal  electrons  are  based  on  expressions  given  by  Schunk  and 
Havs  (1971). 

Figure  10  shows  direct  production  rates  for  three  major  ions  and  the 
heating  rate  resulting  from  the  precipitating  flux  marked  282.3-288.1  s  in 
Figure  11.  This  10  point  average  sample  describes  conditions  near  the  center 
of  the  arc.  The  model  energy  grid  has  19  transporting  energies  from  1800  eV 
to  10  eV  and  7  non-transporting  energies  from  9  eV  to  3  eV.  Electron  in  the 
non-transporting  energy  range  are  assumed  to  deposit  their  energy  locally. 
Fluxes  at  energies  greater  than  1800  eV  are  continued  with  a  power  law 
dependence  to  5  keV.  Electrons  at  the  high  altitude  boundary  are  assumed  to 
be  isotropic  over  the  upper  hemisphere;  zero  upward  flux  through  the  lower 
boundary  at  105  km  altitude  completes  the  boundary  conditions.  Curves  (a) 
and  (b)  in  Figure  10  refer  to  normalized  heating  rates  corresponding  to  the 
two  ambient  electron  density  profiles  (constant  in  time,  cf.  Figure  12). 
Differences  in  production  rates  for  the  assumed  profiles  are  less  than  1%. 
Energy  conservation  tests,  which  compare  the  net  energy  into  the  column 
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with  that  expended  in  ionization  and  excitation,  agree  within  6%  for  case  (a) 
and  2%  for  case  (b). 

Chemistry  model  results  over  a  5  min  interval  using  production  rates 
calculated  for  the  flux  inside  the  arc  are  shown  in  Figure  12.  The  initial 
electron  density  profile  at  t=0  is  a  model  calculation  with  a  representative  flux 
near  the  arc  shown  in  Figure  11.  Also  indicated  in  Figure  12  are  the  ambient 
electron  density  profiles  referred  to  in  Figure  10.  Although  the  model  is 
internally  self-consistent  the  normalized  heating  rate  and  production  rates 
are  held  constant  during  the  calculation.  Since  Figure  10  indicates  that  the 

heating  rate  depends  significantly  on  the  electron  density,  the  corresponding 
Te  profiles  are  markedly  different;  the  electron  density  profiles,  however,  are 
more  nearly  alike  (2-4%  in  the  F-region,  6-10%  in  the  E-region)  and  are  nor 
shown  in  the  figure  for  clarity.  Electric  conductivities  (Figure  13)  are 
determined  from  expressions  given  in  Vallance-Jones  (1974).  Height 
integrated  Pedersen  conductivity  at  5  min  range  from  .76  to  .80  (mho  m‘l) 
using  curves  (b)  and  (a),  respectively. 

Electron  densities  after  2-5  min  in  Figure  12  approximate  the  measured 
values  at  400  km  altitude;  model  Te,  however,  is  600-1000  K  less  than  measured. 
Although  F-region  densities  increase  rapidly  there  is  little  change  at  lower 
altitudes.  This  is  a  response  to  the  large  increase  in  flux  at  energies  less  than  1 
keV;  fluxes  at  energies  greater  than  1  keV  are  very  similar  both  inside  and 
outside  the  arc.  Even  if  the  15  s  profile  were  taken  as  the  initial  condition, 
when  F-region  densities  were  nearly  ten  times  greater  than  at  t=0,  the  final 
profile  after  several  minutes  have  passed  would  be  nearly  identical. 

Although  the  model  electron  density  compares  reasonably  well  with  the 
data,  the  electron  temperature  agreements  is  less  satisfactory.  Because  the 
only  source  now  incorporated  in  the  model  is  collisional  heating  by  energetic 
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electrons,  this  discrepancy  may  in  part  stem  from  other  types  of  heating. 

More  fundamentally,  heat  conduction  is  only  part  of  the  full  energy  balance 
equation  tSehunk  auJ  Walker.  1970);  terms  representing  other  forms  of  heat 
transport  could  modify  the  results.  For  example,  Schunk  et  al.  (1986) 
demonstrate  that  thermoelectric  cooling  by  field-aligned  currents  becomes 
significant  only  for  rather  large  return  currents  of  thermal  electrons  on  the 
order  of  10  flA  m'2.  Weber  et  al.  (1988)  estimate  maximum  field-aligned 
currents  of  about  5  |lA  m'2  within  the  arc,  suggesting  a  minimal  contribution 
from  this  effect. 

Densities  associated  with  an  arc  lasting  a  few  minutes  to  several  tens  of 
minutes  depend  both  on  the  initial  ionosphere  and  the  duration  of  enhanced 
production.  Although  direct  production  by  precipitating  electrons  leads  to  a 
rapid  buildup  of  ionization,  equilibrium  resulting  from  chemistry  and 
diffusion  is  reached  slowly  (several  tens  of  minutes  in  the  F-region  vs. 
minutes  in  the  E-region).  Starting  with  rather  artificial  initial  conditions,  the 
next  two  examples  (Figures  14  and  15)  illustrate  how  the  electron  density 
profile  develops  when  a  flux  like  that  observed  during  PIIE  is  turned  on 
continuously  for  30  minutes. 

In  Figure  14,  0+  ions  alone  are  present  initially  with  a  constant  height 
profile  of  10^  cm'3.  From  120  to  200  km  altitude  production  and  loss  quickly 
balance  and  equilibrium  densities  are  reached  after  10  minutes  or  so.  Below 
120  km  production  is  slight  and  the  density  actually  decreases.  Ion  production 
dominates  above  200  km,  causing  a  steady  increase  in  electron  density  that  is 
limited  by  0+  diffusion  until  equilibrium  is  reached  several  hours  later. 

Starting  with  an  0+  density  of  10^  cm' 3  ,  Figure  15  indicates  that  the  final 
density  profile  after  30  minutes  is  not  different  than  seen  in  Figure  14,  except 
at  the  lowest  altitudes.  From  130  to  200  km  the  equilibrium  densities  are  nearly 
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equal;  above  200  km  the  profiles  are  similar  but  lag  behind  those  in  Figure  14 
because  of  the  smaller  initial  density.  However,  on  time  scales  longer  than  5 
minutes,  tne  resulting  densities  above  iOO  km  agree  tu  within  20%  even  when 
drastically  different  initial  conditions  arc  assumed. 

Figure  16  shows  what  happens  if  the  flux  is  turned  on  for  15  minutes  and 
then  suddenly  switched  off,  after  which  direct  production  ceases  and 
ionization  decays,  subject  only  to  local  chemistry  and  diffusion.  Following  an 
abrupt  drop  caused  by  rapid  losses,  the  electron  density  slowly  decays, 
especially  at  and  above  the  altitude  peak.  Diffusion  flattens  the  profile  and 
maintains  the  density  above  400  km  where  losses  by  recombination  are 
minimal. 

For  comparison.  Figures  17  and  18  show  electron  density  profiles  generated 
by  fluxes  measured  outside  the  main  auroral  arc.  Initial  0+  densities  are  104 
and  10^  cm' ^  as  before.  Here  the  initial  density  conditions  clearly  affect  the 
profiles  even  after  several  minutes  have  elapsed.  Only  in  the  E-region  where 
local  chemistry  is  relatively  fast  do  the  densities  approach  the  same 
equilibrium,  and  neither  case  approximates  the  densities  observed  in  the 
downleg  segment.  Since  the  electron  energy  flux  measured  from  420  km  to  290 
km  nowhere  approaches  the  more  energetic  flux  seen  earlier,  it  seems  likely 
that  the  observed  electron  density  results  from  a  combination  of  ionization 
remaining  after  local  production  ceased  and  that  brought  in  by  convection. 

A  future  goal  of  this  study  is  to  find  the  relationship  between  large  scale 
features  in  the  polar  cap  ionosphere  and  scintillations  of  satellite  to  ground 
radio  transmissions.  Scintillations  are  caused  by  plasma  density  irregularities 
with  scale  sizes  from  100  m  to  10  km  (Weber  and  Buchau.  1985).  These 
irregularities  have  been  associated  with  plasma  instabilities  generated  by 
density  gradients  in  enhanced  plasma  regions  (blobs)  hundreds  of  kilometers 


across  or  with  velocity  shears  (Basu  et  al..  1986).  Hence,  density  features 
resulting  from  direct  production  or  transport  play  an  important  role  in  the 
problem.  Tor  example,  observed  changes  in  scintillation  behavior  with  the 
solar  cycle  (Aarons.  1982)  may  be  connected  with  the  way  these  features  are 
generated.  A  model  can  probe  how  similar  fluxes  affect  neutral  atmospheres 
with  different  exospheric  temperatures.  In  Figure  19,  a  cold  atmosphere  like 
that  found  during  solar  minimum  produces  a  slightly  higher  density  at  lower 
altitudes.  A  hot  neutral  atmosphere  has  peak  density  above  300  km  and  a  more 
gradual  decrease  toward  higher  altitudes  (Figure  20).  Differences  between 
profiles  may  help  identify  which  processes  generate  irregularities  and  the 
scale  sizes  likely  to  be  observed. 
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Table  1 


Chemical  reactions  used  in  the  ionospheric  chemistry  model  (taken  from 
the  model  input  file).  Reactions  are  grouped  by  production  and  loss  for  each 
species.  Rates  equal  to  -1.00E  00  indicate  Te  dependence  interpolated  from  a 
separate  table.  Reactions  with  asterisks  indicate  additions  or  modifications  to 
the  original  set. 

Rate  Reaction 


NO+ 

10  01 
PROD 


13  03- 

•1.00E  00 

02 

1.00 

N2  + 

+ 

0 

- > 

NO+ 

+ 

N 

12  02 

2.00E-10 

00 

N+ 

+ 

02 

- > 

NO+ 

+ 

N 

06  01 

5.00E-16 

00 

02  + 

+ 

N2 

- > 

NO+ 

+ 

NO 

06  08 

4.40E-10 

00 

02  + 

+ 

NO 

- > 

NO+ 

+ 

02 

06  09 

1.80E-10 

00 

02  + 

+ 

N 

- > 

NO+ 

+ 

0 

07  01- 

■1.00E  00 

01 

1.00 

0+ 

+ 

N2 

- > 

NO+ 

+ 

N 

12  08 

8.00E-10 

00 

N+ 

+ 

NO 

- > 

NO+ 

+ 

N 

13  08 

3.30E-10 

00 

N2  + 

+ 

NO 

- > 

NO+ 

+ 

N2 

07  08- 

■l.OOE  00 

18 

1.00 

0+ 

+ 

NO 

- > 

NO+ 

+ 

0 

14  08 
LOSS 

1.20E-C9 

00 

0+( 2D) 

+ 

NO 

- > 

NO+ 

+ 

0 

0.00E 

00  SPONTANEOUS  DECAY  RATE 

05  04- 
02  + 

05  05 
PROD 

•1.00E  00 

03 

NO+ 

+ 

EL 

- > 

N 

+ 

0 

12  02 

4.00E-10 

00 

N+ 

+ 

02 

- > 

02  + 

+ 

N 

07  02- 

■l.OOE  00 

05 

1.00 

0+ 

+ 

02 

- > 

02  + 

+ 

0 

13  02- 

■1.00E  00 

04 

1.00 

N2  + 

+ 

02 

- > 

02+ 

+ 

N2 

14  02 

1.00E-10 

00 

0+ ( 2D) 

+ 

02 

- > 

02  + 

+ 

0 

15  02 

3.00E-10 

00 

0+( 2P) 

+ 

02 

- > 

02  + 

+ 

0 

LOSS 

0.00E 

00 

06  01 

5.00E-16 

00 

02  + 

+ 

N2 

- > 

NO+ 

+ 

NO 

06  08 

4. 40E-10 

00 

02  + 

+ 

NO 

- > 

NO+ 

+ 

02 

06  09 

1.80E-10 

00 

02  + 

+ 

■  N 

- >  NO+ 

+  c 

) 

06  10 

4.00E-10 

00 

02  + 

+ 

N(  2D) 

- > 

N+ 

+ 

02 

06  04- 

•l.OOE  00 

06 

02  + 

+ 

EL 

- > 

0 

+ 

0 

4  5 


05  04-1. 00E  00 

03 

0.20 

N0+ 

+ 

EL 

- > 

N 

+ 

0 

12  03  2.20E-12 

00 

N+ 

+ 

0 

- > 

N 

+ 

0+ 

12  08  8.00E-10 

00 

N+ 

+ 

NO 

- > 

N 

+ 

N0+ 

10  01  1.60E-14 

00 

N(  2D) 

+ 

N2 

- > 

N 

+ 

N2 

10  03  2.00E-13 

00 

N(  2D) 

+ 

0 

- > 

N 

+ 

0 

10  04-1. OOE  00 

12 

1.00 

N {  2D) 

+ 

EL 

- > 

N 

+ 

EL 

19  03  5.00E-11 
LOSS 

O.OOE  00 

00 

N2 ( A3S ) 

+ 

0 

- > 

N 

+ 

NO 

09  02-1. OOE  00 

11 

N 

+ 

02 

- > 

NO 

+ 

0 

09  08-1. OOE  00 

10 

N 

+ 

NO 

- > 

N2 

+ 

0 

09  06  1.20E-10 
0+ 

07  05 

PROD 

00 

N 

+ 

02  + 

- > 

N0+ 

+ 

0 

14  03  O.OOE-11 

00 

0: 

(2D) 

+ 

0 

- > 

0+ 

+ 

0 

14  04-1. OOE  00 

09 

0.50 

0+ { 2D) 

+ 

EL 

- > 

0+ 

+ 

EL 

15  04-1. OOE  00 

09 

0.50 

0+ ( 2P ) 

+ 

EL 

- > 

0+ 

+ 

EL 

13  03-1. OOE  00 

07 

1.00 

N2  + 

+ 

0 

- > 

0+ 

+ 

N2 

12  03  2.20E-12 

00 

N+ 

+ 

0 

- > 

0+ 

+ 

N 

12  02  2.00E-11 

00 

N+ 

+ 

02 

- > 

0+ 

+ 

NO 

14  01  2.00E-11 
LOSS 

O.OOE  00 

00 

0+ ( 2D) 

+ 

N2 

- > 

0+ 

+ 

N2 

07  01-1. OOE  00 

01 

0+ 

+ 

N2 

- > 

N0+ 

+ 

N 

07  02-1. OOE  00 

05 

0+ 

+ 

02 

- > 

02  + 

+ 

0 

07  08-1. OOE  00 

18 

0+ 

+ 

NO 

- > 

N+ 

+ 

02 

07  10  1.30E-10 

00 

0+ 

+ 

N(  2D) 

- > 

N+ 

+ 

0 

07  04  1.60E-12 
NO 

00  01 

LOSS 

O.OOE  00 

08  06  O.OOE  QO 
N(4S) 

11  03 

PROD 

00 

00 

0+ 

+ 

EL 

- > 

0 

+ 

( 7774A) 

13  04-1. OOE  00 

08 

1.00 

N2+ 

+ 

EL 

- > 

N 

+ 

N(  2D) 

13  03-1. OOE  00 

02 

0.05 

N2  + 

+ 

0 

- > 

N 

+ 

N0+ 

12  02  3.00E-10 

00 

N+ 

+ 

02 

- > 

N 

+ 

02  + 

07  01-1. OOE  00 

01 

1.00 

0+ 

+ 

N2 

- > 

N 

+ 

N0+ 

46 


N(  2D) 
05  07 
PROD 


13  04-1. 00E  00 

08 

1.00 

N2  + 

+ 

EL 

- > 

N(  2D) 

+ 

N 

13  03-1. 00E  00 

02 

0.95 

N2  + 

+ 

0 

- > 

N(  2D) 

+ 

N0+ 

12  02  3.90E-10 

00 

N+ 

+ 

02 

- > 

N(  2D) 

+ 

02  + 

05  04-1. OOE  CO 

03 

0.80 

N0+ 

+ 

EL 

- > 

N(  2D ) 

+ 

0 

19  03  5.00E-11 
LOSS 

1.10E-05 

00 

N2 ( A3S ) 

+ 

0 

- > 

N(  2D) 

+ 

NO 

10  01  1.60E-14 

00 

N (  2D) 

+ 

N2 

- > 

N 

+ 

N2 

10  02-1. OOE  00 

19 

N(  2D) 

+ 

02 

- > 

N 

+ 

02 

10  08  7.00E-11 

00 

N{  2D) 

+ 

NO 

- > 

N 

+ 

NO 

10  06  4.00E-10 

00 

N(  2D) 

+ 

02  + 

- > 

N+ 

+ 

02 

10  03  4.00E-13 

00 

N(  2D) 

0 

- > 

N 

+ 

0 

10  04-1. OOE  00 

12 

N{  2D) 

+ 

EL 

- > 

N 

+ 

EL 

10  07  1.30E-10 
02 ( 1DELG ) 

00  03 

LOSS 

3.70E-04 

00 

N(  2D) 

+ 

0+ 

- > 

N+ 

+ 

0 

11  09-1. OOE  00 

13 

02 (IDG) 

+ 

N 

- > 

02 

+ 

N 

11  03  1.00E-17 

00 

02(1DG) 

+ 

0 

- > 

02 

+ 

0 

11  02-1. OOE  00 
N+ 

02  05 

PROD 

14 

02  (IDG) 

+ 

02 

- > 

02 

+ 

02 

06  10  4.00E-10 

00 

02  + 

+ 

N (  2D) 

- > 

N+ 

+ 

02 

07  10  1.30E-10 
LOSS 

0.00E  00 

00 

0+ 

+ 

N(  2D) 

- > 

N+ 

+ 

0 

12  02  2.00E-11 

00 

N+ 

+ 

02 

- > 

0+ 

+ 

NO 

12  02  2.00E-10 

00 

N+ 

+ 

02 

- > 

N0+ 

+ 

0 

12  02  4.00E-10 

00 

N+ 

+ 

02 

- > 

02  + 

+ 

N 

12  03  2.20E-12 

00 

N+ 

+ 

0 

- > 

0+ 

+ 

N 

12  08  8.00E-10 
N2  + 

03  05 

PROD 

00 

N+ 

+ 

NO 

- > 

N0+ 

+ 

N 

16  01  2.50E-10 
14  01  1.00E-10 

00 

00 

02+ ( A4P ) 
0+( 2D) 

+ 

N2 

+  N2 

- >  N2  + 

- >  N2+  + 

+ 

0 

02 

15  01  4.80E-10 
LOSS 

0.00E  00 

00 

0+( 2P) 

+ 

N2 

- > 

N2  + 

+ 

0 

13  02-1. OOE  00 

04 

N2  + 

+ 

02 

- > 

02+ 

+ 

N2 

13  03-1. OOE  00 

02 

N2  + 

+ 

0 

- > 

N0+ 

+ 

N(  2D) 

13  08  3.30E-10 

00 

N2  + 

+ 

NO 

- > 

N0+ 

+ 

N2 

13  04-1. OOE  00 

08 

N2  + 

+  E 

- ; 

>  N 

f  N 

13  03-1. OOE  00 

07 

1.00 

N2+  + 

0 

- >  0+ 

+N2 

0+( 2D) 

00  05 

LOSS 

7.70E-05 

14  01  1.00E-1C 

00 

0+ ( 2D)  + 

N2 

- > 

N2  + 

+ 

0 

14  02  1.00E-10 

00 

0+ ( 2D)  + 

02 

- > 

02  + 

+ 

0 

14  03  0.00E-11 

00 

0+ ( 2D)  + 

0 

- > 

0+ 

+ 

0 

14  04-1. 00E  00 

09 

0+ ( 2D )  + 

EL 

- > 

0+ 

+ 

EL 

14  08  1.20E-09 

00 

0+ ( 2D)  + 

NO 

- > 

N0+ 

+ 

0 

0+ ( 2P ) 

00  03 

LOSS 

2.00E-01 

15  01  4.80E-10 
15  02-1. OOE  00 

00 

05 

0+ ( 2P ) 

0+ ( 2P )  + 

+  N2 

02 

- >  N2+  + 

— >  02  + 

0 

+ 

0 

15  04-1. OOE  00 

09 

0+ ( 2P )  + 

EL 

- > 

0+ 

+ 

EL 

02+ ( A4PI ) 

00  04 

LOSS 

O.OOE  00 

16  01  2.50E-10 

00 

02+ ( A4P )  + 

N2 

- > 

N2  + 

+ 

02 

16  02  3.00E-10 

00 

02+ ( A4P )  + 

02 

- > 

02  + 

+ 

02 

16  03  1.00E-10 

00 

02+ ( A4P )  + 

0 

- > 

0+ 

+ 

02 

16  04  1.00E-07 

00 

02+(A4P)  + 

EL 

- > 

0 

+ 

0 

0(  ID) 

05  03 

PROD 

10  08  2.00E-11 

00 

N( 2D)  + 

NO 

- > 

0(  ID) 

+ 

N2 

10  03  2.00E-13 

00 

N( 2D)  + 

0 

- > 

0  ( ID ) 

+ 

N 

10  02  4.80E-12 

00 

N  (  2D )  + 

02 

- > 

0(  ID) 

+ 

NO 

06  04-1. OOE  00 

06 

1.00 

02+  + 

E 

- > 

0(  ID) 

+ 

0 

18  03-1. OOE  00 

20 

2.00 

0(1S)  + 

0 

- > 

0  ( ID ) 

+ 

0  ( ID ) 

LOSS 

6.74E-03 

17  01  2.30E-11 

00 

0( ID)  + 

N2 

- > 

0 

+ 

N2 

17  08  3.50E-11 

00 

0( ID)  + 

NO 

- > 

0 

+ 

NO 

17  02  4.10E-11 

00 

0( ID)  + 

02 

- > 

0 

+ 

02 

0(  IS) 

04  04 

PROD 

12  02  1.90E-10 

00 

N+  + 

02 

- > 

N0+ 

+ 

0 

06  04-1. OOE  00 

06 

0.08 

02+  + 

E 

- > 

0(1S) 

+ 

0 

19  03  3.60E-11 

00 

N2 ( A3S ) (1)+ 

0 

- > 

0(1S) 

+ 

N2 

06  09  2.50E-11 

00 

02+  + 

N(4S) 

- > 

0(1S) 

+ 

NO 

LOSS 

1.41E  00 

18  03  2.00E-14 

00 

0(1S)  + 

0 

- > 

0  ( ID ) 

+ 

0(  ID) 

18  08-1. OOE  00 

15 

0(1S)  + 

NO 

- > 

0 

+ 

NO 

18  02-1. OOE  00 

16 

0 ( IS )  + 

02 

- > 

0 

+ 

02 

18  11  1.70E-10 

00 

0(1S)  + 

02  (IDG) 

- > 

0 

+ 

02 

48 


N2 ( A3SIG) 
00  04 
LOSS 
5.30E-01 


19  02  2.00E-11 

00 

N2 ( A3S ) 

+ 

02 

- >  N2 

+ 

02 

19  03  3.40E-11 

00 

N2 ( A3S ) ( 1 )  + 

0 

- >  N2 

+ 

0 

19  08  4.00E-11 

00 

N2 ( A3S ) ( 1 ) + 

NO 

- >  N2 

+ 

NO 

19  09  5.00E-11 
02 ( B1SIG ) 

00 

N2 ( A3S ) 

+ 

N 

- >  N2 

+ 

N 

01  02 

PROD 

17  02  5.00E-11 
LOSS 

8.30E-02 

00 

0(  ID) 

+ 

02 

- >  0 

+ 

02 ( BIS ) 

20  02-1. 00E  00 

17 

02  ( BIS ) 

+ 

02 

- >  02 

+ 

02 

20  01  1.00E-17 

00 

02 ( BIS ) 

+ 

N2 

- >  02 

+ 

N2 

4  9 


Table  2 


Comparison  of  model  Te  with  analytic  expression  (both  at  steady  state). 
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initial  Te 

(K)  A 

S 

model  Te 

(K)  analytic  1 

600 

2048 

1.46(14) 

7.20(0) 

1494 

1495 

500 

2048 

1.46(14) 

5.10(1) 

1488 

1489 
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2033 

1.44(14) 

1.02(2) 
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1475 
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1.38(14) 
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1.28(14) 
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6.54(1) 

1055 

1056 

200 

1360 
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-2.46(2) 
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1.84(13) 

-5.12(2) 
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362 

50 


Table  3 


Tabular  Production  rates  (cm"3  s"1)  for  PIIE  flux  (full  model  results  in 


parentheses 

when 

different). 
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